
%Data = csvread('Rawdata.csv',1,0);
Data = readtable('Rawdata_1sector.csv');

data.ind = 1;
data.NIC = 99;
data.N_i = Data.N_i;
data.N_f = Data.N_f;
data.X_i = Data.X_i*1e5;
data.X_f = Data.X_f*1e5;
data.L_i = Data.L_i*1e3;
data.L_f = Data.L_f*1e3;
data.L_c = Data.L_c*1e3;
%data.l_i = Data.LnL_i;
%data.l_f = Data.LnL_f;
data.Vl_i = Data.VL_i;
data.Vl_f = Data.VL_f;
data.Lbar_i = Data.Lbar_i;
data.Lbar_f = Data.Lbar_f;
%data.lbar_i = Data.LnLbar_i;
%data.lbar_f = Data.LnLbar_f;

pm.kappa = (data.X_i + data.X_f)./sum((data.X_i + data.X_f));
pm.t = Data.Tax./100;

% Two ways to measure the ratio of contract workers
pm.ratio = Data.Frac_contract./(1-Data.Frac_contract); %average across firms
%pm.ratio = Data.Frac_contract2./(1-Data.Frac_contract2); %average across firms
%data.ratio = data.L_c./(data.L_f-data.L_c); %fraction of workers on a contract

%Calibrated parameters of the model
data.lbar_i = log(Data.Lbar_i);
data.lbar_f = log(Data.Lbar_f);

pm.L= sum(data.L_i + data.L_f);
pm.I = size(Data,1); %total number of sectors
pm.nu = 3;
pm.rho = 0.738;
%pm.phi = 1;s
data.coeff_2 = -0.0100236;

%{
Source:
https://censusindia.gov.in/Tables_Published/A-Series/A-Series_links/t_00_009.aspx

Total workers in AP are around 35 million. Total non-agricultural workers
are 11.5 million. For now, we can take those to be the "potential entrants"
for the number of firms
%}
pm.M = repmat(11.5*1e6 - sum(data.N_i+data.N_f) - pm.L,pm.I,1);
%pm.M = repmat(5.*max((data.N_i+data.N_f)),pm.I,1);
%pm.M = 10.*(data.N_i+data.N_f);

pm.mintolp = 0.005;
%B = (1 + data.ratio.^((1-pm.nu)./pm.nu)).^(1./(1-pm.nu));


%% ESTIMATION PROCEDURE 

fprintf('Estimating parameters \n');
Param= readmatrix('Sigmastar_1sector.csv');

%Data Moments
M_data = [mean(Data.VL_i);mean(Data.VL_f); Data.Lbar_i;Data.Lbar_f; data.coeff_2];

%STEP 1: Calculate the matrix without any weighting parameter
L_unwt = @(param) (Algorithm(param,pm,Data,1)./M_data-1)' * (Algorithm(param,pm,Data,1)./M_data-1);

P0 = [  Param(3);Param(4);...
        Param(5);...
        Param(7); ...
        Param(6)];
%P0 = [  0.2;0.2; 1.1; 0.8; 1];

lb =    [0.1;0.1;...
        ones(pm.I,1);...
        zeros(pm.I,1); ...
        0.8];
ub = [  5;5;...
        10.*ones(pm.I,1);...
        10.*ones(pm.I,1); ...
        1.2];
    
sigoptions = optimoptions('fmincon','Display','iter','UseParallel',false, 'MaxIter',5500000,'OptimalityTolerance',1e-4);
%sigoptions = optimoptions('fmincon','Display','iter', 'MaxIter',5500000,'MaxFunEvals',5500000);
paramstar= fmincon(L_unwt, P0,[],[],[],[],lb,ub,[],sigoptions);

pm.sigma = paramstar(1);
pm.se = paramstar(2);       
pm.theta = paramstar(3);
pm.bc = paramstar(4);
pm.phi = paramstar(5);

zr =logninv(1-Data.N_f./pm.M,0,pm.sigma);
pm.bp = zeros(size(pm.bc));
for r = 1:pm.I
    xbar = integral(@(x) x.^pm.phi.*lognpdf(x,0,pm.sigma),zr(r),Inf) ...
                ./(1-logncdf(zr(r),0,pm.sigma));
    Ee= integral(@(e) e.^pm.phi.*lognpdf(e,0,pm.se),-Inf,Inf);
    pm.bp(r) = pm.ratio.^(1./pm.nu).*pm.bc./(xbar.*Ee);
end

%Model Fit using the data
ze =logninv(1-(Data.N_f+Data.N_i)./pm.M,0,pm.sigma);
zr = logninv(1-Data.N_f./pm.M,0,pm.sigma);
[p,w,pm.E_i,pm.E_r]= PriceWages(pm,ze,zr);

writematrix([data.ind data.NIC ...
            pm.sigma pm.se ...
            pm.theta pm.phi pm.bc pm.bp ...
            pm.E_i pm.E_r ...
            p w Data.N_f+Data.N_i],'Sigmastar_1sector.csv','delimiter',',');
fprintf('\n \n Done: Estimation \n \n');


%% MODEL FIT
Param = readmatrix('Sigmastar_1sector.csv');

pm.sigma = Param(3);
pm.se = Param(4);
pm.theta = Param(5);
pm.phi = Param(6);
pm.bc = Param(7);
pm.bp = Param(8);
pm.E_i = Param(9);
pm.E_r = Param(10);
pm.price = Param(11);
pm.wage = Param(12);

%Model Fit and Baseline scenario: find the baseline wages in the economy before we change
%the bc/bp ratios.
i = 1;
adjfactor = Inf;
[p,w,ZMat,ZstarMat,LMat,YMat,NMat,Income,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
Sim.VL_i = VMat(1);
Sim.VL_f = VMat(2);

ze = ZstarMat(:,1);
zr = ZstarMat(:,2);

Mlbar = Sim_Firmsize(pm,p,w,ze,zr);
Sim.Lbar_i = Mlbar(1:pm.I);
Sim.ratio = zeros(pm.I,1);
for r = 1:pm.I
    xbar = integral(@(x) x.^pm.phi.*lognpdf(x,0,pm.sigma),zr(r),Inf) ...
                ./(1-logncdf(zr(r),0,pm.sigma));
    Ee= integral(@(e) e.^pm.phi.*lognpdf(e,0,pm.se),-Inf,Inf);
    Sim.ratio(r) = (pm.bc./(xbar.*Ee.*pm.bp)).^(-pm.nu);
end

zgrid = linspace(1,2,10)'.*zr;
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
Sim.coeff = polyfit(lnL,cr,2);

%Simulating the fraction of contract workers by firm size 
gpts = 1000;
x = repmat(linspace(0.1,5.2,gpts),pm.I,1);
e = lognrnd(0,pm.se,gpts);

L_f = zeros(pm.I,gpts);
L_c = zeros(pm.I,gpts);
L_p = zeros(pm.I,gpts);

for r = 1:pm.I
    
    L_f(r,:) = integral(@(e) l_f(pm,x,e,p(r),w,pm.t(r),r)...
        .*lognpdf(e,0,pm.se),-Inf,Inf,'ArrayValued',true);

    L_c(r,:) = integral(@(e) ...
            (w.*pm.bc(r)./((x.*e).^pm.phi.*wf(pm,x,e,w,r))).^(-pm.nu) .*L_f(r,:) ...
        .*lognpdf(e,0,pm.se),-Inf,Inf,'ArrayValued',true);

    L_p(r,:) = integral(@(e) ...
            (w.*pm.bp(r)./wf(pm,x,e,w,r)).^(-pm.nu) .*L_f(r,:) ...
        .*lognpdf(e,0,pm.se),-Inf,Inf,'ArrayValued',true);
end

Sim.Lcorr = [L_c./(L_c+L_p); L_c + L_p ]';
writematrix([L_c./(L_c+L_p); L_c + L_p ]','Lcorr_1sector.csv','delimiter',',');

Sim.Lbar_f = mean(L_f(L_f<1000));

%Outsheeting Model Fit
writematrix([data.ind data.NIC  ...
             NMat(:,1) Data.N_i ...
             NMat(:,2) Data.N_f ...
             LMat(:,1) data.L_i ...
             LMat(:,2) data.L_f ...
             Sim.Lbar_i Data.Lbar_i ...
             Sim.Lbar_f Data.Lbar_f  ...
             Sim.ratio./(1+Sim.ratio) pm.ratio./(1+pm.ratio) ...
             repmat([Sim.VL_i mean(Data.VL_i) ...
             Sim.VL_f mean(Data.VL_f)],pm.I,1) ...
             Sim.coeff(1) data.coeff_2],'Modelfit_1sector.csv','delimiter',',');


%writematrix([data.nic p LMat YMat NMat],'Modelfit.csv','delimiter',',');
%writematrix([data.nic pm.kappa pm.bc pm.bp p repmat([w i Income],pm.I,1) ZMat ZstarMat LMat YMat NMat],['Counterfactuals_' num2str(i) '.csv'],'delimiter',',');
base.bc = pm.bc;
base.bp = pm.bp;
pm.wbase = w;

fprintf('\n \n Done: Model Fit \n \n');


%% COUNTERFACTUALS


igrid = [0 19];

%Change in bc x
for inum = 1:size(igrid,2)
    i = igrid(inum);
    bc = (1+0.01.*i).*base.bc;
    
    %presence of informal sector
    [p,w,ZMat,ZstarMat, LMat,YMat,NMat,Income,VMat,TFP,Cratio] = Counterfactuals(pm,bc,base.bp);

    writematrix([data.NIC pm.kappa bc base.bp p Cratio repmat([w i Income TFP],pm.I,1)  ZMat ZstarMat LMat YMat NMat],...
                ['Counterfactuals_1s_bc_' num2str(inum) 'A.csv'],'delimiter',',');
    fprintf('Done with %d \n',inum);
end


%Change in bp 
for inum = 2:size(igrid,2)
    i = igrid(inum);
    bp = (1-0.0065.*i).*base.bp;

    %presence of informal sector
    [p,w,ZMat,ZstarMat, LMat,YMat,NMat,Income,VMat,TFP,Cratio] = Counterfactuals(pm,base.bc,bp);
    writematrix([data.NIC pm.kappa base.bc bp p Cratio repmat([w i Income TFP],pm.I,1)  ZMat ZstarMat LMat YMat NMat],...
                ['Counterfactuals_1s_bp_' num2str(inum) 'A.csv'],'delimiter',',');
end

fprintf('\n \n Done: Counterfactuals \n \n');
